#include "VarianceGammaDistribution.h"
#include <math.h>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <gsl/gsl_specfunc.h>
#include <string>

VarianceGammaDistribution::VarianceGammaDistribution(double mu, double lambda, double alpha, double beta) 
	: _mu(mu), _lambda(lambda), _alpha(alpha), _beta(beta), _gamma(sqrt(alpha*alpha - beta*beta)), _pi(boost::math::constants::pi<double>()),
		_quantileBounds(mu, lambda, alpha, beta), _asympF(mu, lambda, alpha, beta, 20) 
{
	using namespace boost::math;

	if(!validParameters(mu, lambda, alpha, beta))
		throw std::exception("VG parameters invalid");

	_mean = _mu + ((2 * _beta * _lambda)/(_gamma * _gamma));
	_variance = ((2 * _lambda * (1+(2 * _beta * _beta)/(_gamma * _gamma)))/(_gamma * _gamma));
	_normalizingConstant =  pow(_gamma, 2 * _lambda)/(pow(2 * _alpha, _lambda - 0.5) * sqrt(_pi) * tgamma(_lambda));
	_support = findSupport();
}

VarianceGammaDistribution::~VarianceGammaDistribution(void)
{
}

double VarianceGammaDistribution::density(double x) const
{
	using namespace boost::math;
	
	if(x == _mu)
		return std::numeric_limits<double>::infinity();

	/*gsl_sf_bessel_Knu(_lambda-0.5, _alpha*abs(x - _mu))*/
	double result;
	std::string error;
	
	try{
		result = _normalizingConstant * pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu));
	}catch(...){
		std::stringstream errMsg;
		errMsg << "failed to evaluate density at x = " << x;
		throw std::runtime_error(errMsg.str());
	}
	/*
	try
	{
		result = _normalizingConstant * pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu));

	}catch( std::domain_error e)
	{
		error = e.what();
		result = 0;
	}catch( std::overflow_error e)
	{
		error = e.what();
		result = 0;
	}catch( std::underflow_error e)
	{
		error = e.what();
		result = 0;
	}catch( boost::math::evaluation_error e)
	{
		error = e.what();
		result = 0;
	}
	*/
	return result;
}

std::pair<double,double> VarianceGammaDistribution::quantileBounds( double u ) const
{
	return _quantileBounds(u);
}

std::pair<double, double> VarianceGammaDistribution::findSupport() const
{
	// find 2 points for which the density is greater than eps. 
	double eps = std::numeric_limits<float>::epsilon();
	double stdDev = sqrt(variance());
	double step = stdDev/2;
	double xLeft = _mu - 4*stdDev;
	double xRight = _mu + 4*stdDev;

	const size_t MAX_ITERATIONS = 100000;

	//unsigned int k = 1;
	bool left = false;
	bool right = false;

	// Find break points
	for(size_t n = 1; n < MAX_ITERATIONS; n++)
	{
		if(!left && density(xLeft) <= eps)
			left = true;
		

		if(!right && density(xRight) <= eps)
			right = true;
		
		if(left && right)
			break;

		xLeft = _mu - (n*step);
		xRight = _mu + (n*step);
	}

	if(!(left&&right))
		throw std::logic_error("Could not find the numerical support of the distribution within iteration limit.");

	return make_pair(xLeft, xRight);
}

double VarianceGammaDistribution::cdf( double x )
{
	double xL = getSupport().first;
	double xR = getSupport().second;

	// Integrate density, for the tail regions we will use the 
	// asymptotic expansion of F. 

	if( x <= xL ) 
		return _asympF.LeftTail(x);

	if( xL < x && x <= xR )
		return _asympF.LeftTail(xL) + probability(xL, x);
	
    // i.e x > bR
	return _asympF.RightTail(x);

}

double VarianceGammaDistribution::probability( double a, double b )
{
	//TODO: Add some error handling  
	double result, error;      
	gsl_function F;
	F.function = &ProbabilityDistribution::densityWrapper;
	F.params = this;

	double epsabs = 0;
	double epsrel = 1e-13;

	try
	{
		if(b > _mu){
			double pts[3] = {a, _mu, b};
			gsl_integration_qagp (&F, pts, 3, epsabs, epsrel, 1000, _integrationWorkspace, &result, &error);
		}else
		{
			gsl_integration_qag (&F, a, b, epsabs, epsrel, 1000, GSL_INTEG_GAUSS61, _integrationWorkspace, &result, &error);
		}
	}
	catch (...)
	{
		std::stringstream errMsg;
		errMsg << "failed to determine probability of the event {a < X < b}, where a = " << a << "and b = " << b;
		throw std::runtime_error(errMsg.str());
	}
	


	//printf("result          = % .18f\n", result);
	//printf ("estimated error = % .18f\n", error);
	//printf ("intervals =  %d\n", w->size);

	return result;
}

bool VarianceGammaDistribution::validParameters( double mu, double lambda, double alpha, double beta )
{
	// TODO: parameter validation
	if(abs(beta)>_alpha)
		return false;

	return true;
}